function histonesPlots_CancervsNormal_Fig4B
% Comparison of histone modification data between a normal (GM06990) and cancer cell line (K562), over chromosome 7
%
% SYNTAX: [] = histonePlots_CancervsNormal_Fig4B()
%
% See also assembleHistoneDSforClustering(), assembleHistoneDS()
%  
%    DASMiner: DAS library and browser for Matlab.
%    Diogo Veiga, Jan 2009.

[xmlS, xml] = executeDASCommand('http://hgwdev-gencode.cse.ucsc.edu/cgi-bin/das','command', 'features', 'DSN', 'hg18','chrom','7','start','115404471' ...
    , 'stop', '117281897', 'featuresList',{'encodeSangerChipH3acK562'});
[avg_score_H3ac_K562] = readHistoneXML(xmlS);

[xmlS, xml] = executeDASCommand('http://hgwdev-gencode.cse.ucsc.edu/cgi-bin/das','command', 'features', 'DSN', 'hg18','chrom','7','start','115404471' ...
    , 'stop', '117281897', 'featuresList',{'encodeSangerChipH3ac'});
[avg_score_H3ac] = readHistoneXML(xmlS);

[xmlS, xml] = executeDASCommand('http://hgwdev-gencode.cse.ucsc.edu/cgi-bin/das','command', 'features', 'DSN', 'hg18','chrom','7','start','115404471' ...
    , 'stop', '117281897', 'featuresList',{'encodeSangerChipH3K4me3K562'});
[avg_score_H3K4me3_K562] = readHistoneXML(xmlS);

[xmlS, xml] = executeDASCommand('http://hgwdev-gencode.cse.ucsc.edu/cgi-bin/das','command', 'features', 'DSN', 'hg18','chrom','7','start','115404471' ...
    , 'stop', '117281897', 'featuresList',{'encodeSangerChipH3K4me3'});
[avg_score_H3K4me3] = readHistoneXML(xmlS);

                         
%Transform intensity into [0,1] interval
avg_score = [avg_score_H3ac avg_score_H3ac_K562 avg_score_H3K4me3 avg_score_H3K4me3_K562];
[avg_score] = mapminmax(avg_score',0,1);
avg_score_H3ac = avg_score(1,:);
avg_score_H3ac_K562 = avg_score(2,:);
avg_score_H3K4me3 = avg_score(3,:);
avg_score_H3K4me3_K562 = avg_score(4,:);

save 'histonePlots_data.mat' avg_score_H3ac avg_score_H3ac_K562 ...
                             avg_score_H3K4me3 avg_score_H3K4me3_K562;
                         
%plottools('on')

hs_cytobands = cytobandread('hs_cytoBand.txt');
chromosomeplot(hs_cytobands, 7, 'Orientation', 2, 'Unit', 1);

subplot(4,1,1);
plot(avg_score_H3ac, '-r', 'DisplayName', 'GM06990 H3ac', 'YDataSource', 'avg_score_H3ac');
title('GM06990 H3ac')
box off

subplot(4,1,2); 
plot(avg_score_H3ac_K562,  '-r', 'DisplayName', 'avg_score_H3ac_K562', 'YDataSource', 'avg_score_H3ac_K562');
title('K562 H3ac')
box off

subplot(4,1,3); 
plot(avg_score_H3K4me3, '-b', 'DisplayName', 'avg_score_H3K4me3', 'YDataSource', 'avg_score_H3K4me3');
title('GM06990 H3K4me3')
box off

subplot(4,1,4); 
plot(avg_score_H3K4me3_K562, '-b', 'DisplayName', 'avg_score_H3K4me3_K562', 'YDataSource', 'avg_score_H3K4me3_K562');
title('K562 H3K4me3')
box off



function [avg_score,absPos,nHits]=readHistoneXML(xmlS)

nFeatures = size(xmlS.GFF{1}.SEGMENT{1}.FEATURE,2);
feat = cell(nFeatures,1);

chrEnd = 117281897;
chrStart = 115404471;

absPos = zeros(chrEnd-chrStart+1,1);
for i=1:chrEnd-chrStart+1
   absPos(i) = i; 
end
avg_score = zeros(chrEnd-chrStart+1,1);
nHits = zeros(chrEnd-chrStart+1,1);

for i=1:nFeatures
    
    %chipLabel = xmlS.GFF{1}.SEGMENT{1}.FEATURE{i}.METHOD{1}.CONTENT;
    start = str2double(xmlS.GFF{1}.SEGMENT{1}.FEATURE{i}.START{1}.CONTENT);
    stop = str2double(xmlS.GFF{1}.SEGMENT{1}.FEATURE{i}.END{1}.CONTENT);
    score = str2double(xmlS.GFF{1}.SEGMENT{1}.FEATURE{i}.SCORE{1}.CONTENT);
    
    if (start-chrStart) < 0
        display('out fragment');
    end
    
    for j=start:stop
        
        if (j > chrEnd) , 
            continue; 
        end
        if (j < chrStart), 
            continue; 
        end 
            
       nHits(j-chrStart+1) = nHits(j-chrStart+1)+1;
       avg_score(j-chrStart+1) = (avg_score(j-chrStart+1) + score)/nHits(j-chrStart+1);
    end
end

